Homework 7.1: Cytokines & Homeostasis

a) Nondimensionalization

\begin{align} \cfrac{\mathrm{d}c}{\mathrm{d}t} &= \beta_0 \cfrac{c(m/k)^2}{1+(m/k)^2} - \gamma_c c - \alpha_c m c \\[0.5em] \cfrac{\mathrm{d}m}{\mathrm{d}t} &= I_0 + \beta_1 c - \alpha_0 \cfrac{c (m/k)^2}{1+(m/k)^2} - \gamma_m m \end{align}

$\hspace{15em}\mathrm{Nondimensionalizing... }\\[0.5em]$ \begin{align} \cfrac{\mathrm{d}\tilde{c}}{\mathrm{d}\tilde{t}} &= \beta_0 \cfrac{\tilde{c}(\tilde{m} m_d /k)^2}{1+(\tilde{m}m_d/k)^2} t_d - \gamma_c \tilde{c} t_d - \alpha_c m_d \tilde{m} \tilde{c} t_d \\[0.5em] \cfrac{\mathrm{d}\tilde{m}}{\mathrm{d}\tilde{t}} &= I_0\cfrac{t_d}{m_d} + \beta_1 c_d \tilde{c} \cfrac{t_d}{m_d} - \alpha_0 \cfrac{c_d \tilde{c} (m_d \tilde{m}/k)^2}{1+(m_d \tilde{m}/k)^2}\cfrac{t_d}{m_d} - \gamma_m \tilde{m} t_d \\[1.5em] \end{align}

$\hspace{15em}\mathrm{Setting... }$ \begin{align} \boxed{m_d = k, \hspace{0.5em} t_d = 1/\gamma_m, \hspace{0.5em} c_d = \cfrac{\gamma_m k}{\beta_1}}\\[0.1em] \end{align} \begin{align} \\[0.5em] \cfrac{\mathrm{d}\tilde{c}}{\mathrm{d}\tilde{t}} &= \cfrac{\beta_0}{\gamma_m} \cfrac{\tilde{c}\tilde{m}^2}{1+\tilde{m}^2} - \cfrac{\gamma_c}{\gamma_m} \tilde{c} - \cfrac{\alpha_c k}{\gamma_m}\tilde{m}\tilde{c} \\[0.5em] \cfrac{\mathrm{d}\tilde{m}}{\mathrm{d}\tilde{t}} &= \cfrac{I_0}{k \gamma_m} + \cfrac{\beta_1}{\gamma_m k} c_d \tilde{c} - \cfrac{a_0}{\gamma_m k}c_d \tilde{c} \cfrac{\tilde{m}^2}{1+\tilde{m}^2} - \tilde{m}\\[1.5em] \end{align}

$\hspace{15em}\mathrm{Setting... }$ \begin{align} \boxed{\tilde{\beta_0} = \beta_0/\gamma_m, \hspace{0.5em} \gamma = \gamma_c/\gamma_m, \hspace{0.5em} \tilde{\alpha_c} = \cfrac{\alpha_c k}{\gamma_m} \\[0.5em] \hspace{2em} \tilde{I_0} = \cfrac{I_0}{k \gamma_m}, \hspace{0.5em} \tilde{\alpha_0} = \alpha_0/\beta_1 } \end{align} \begin{align} \\[0.1em] \cfrac{\mathrm{d}\tilde{c}}{\mathrm{d}\tilde{t}} &= \tilde{\beta_0} \cfrac{\tilde{m}^2 \tilde{c}}{1+\tilde{m}^2} - \gamma \tilde{c} - \tilde{\alpha_c} \tilde{m} \tilde{c} \\[0.5em] \cfrac{\mathrm{d}\tilde{m}}{\mathrm{d}\tilde{t}} &= \tilde{I_0} + \tilde{c} - \tilde{\alpha_0} \tilde{c} \cfrac{\tilde{m}^2}{1+\tilde{m}^2} - \tilde{m}\\[1.5em] \end{align} $\hspace{15em}\mathrm{Rewriting... }$ \begin{align} \boxed{ \cfrac{\mathrm{d}c}{\mathrm{d}t} = {\beta_0} \cfrac{m^2 c}{1+m^2} - \gamma c - \alpha_c m c \\[0.5em] \cfrac{\mathrm{d}m}{\mathrm{d}t} = I_0 + c - \alpha_0 c \cfrac{m^2}{1+m^2} - m } \end{align}

b) Nullclines for $\cfrac{\mathrm{d}c}{\mathrm{d}t} = 0$

  1. \begin{align} \cfrac{\mathrm{d}c}{\mathrm{d}t} = c\left(\beta_0 \cfrac{m^2}{1+m^2}-\gamma - \alpha_c m \right) = 0 \end{align}

No matter what the concentration of $m$, if $c = 0$, then c is not changing. Thus in the $c-m$ plane, this is a vertical line at $c = 0$, and $m$ takes on the entire range of real (positive) numbers. Note that when we plug in $c = 0$ to dm/dt, then we get $m = I_0$.

  1. When $c$ is nonzero, for the expression to hold, the second term must be zero, which means we can solve for a cubic of $m$ to give: \begin{align} \boxed{\alpha_c m^3 + (\gamma - \beta_0) m^2 + \alpha_c m + \gamma = 0} \end{align} By DesCartes's Rule of Signs, zero sign changes means zero positive roots, and two sign changes means a possible maximum of two positive roots. Since $\alpha_c, \gamma$ are positive, only the second term has potential to be negative, i.e. when $\gamma - \beta_0 < 0$, when $\beta_0$ > $\gamma$. This means that the dimensional parameters $\beta_0 / \gamma_m > \gamma_c / \gamma_m$, or when $\beta_0 > \gamma_c$. Physically, this means that the production rate of $c$ responds faster than the baseline degradation of $c$. Since $c$ appears nowhere in this term, that means when we find values for $m$ that send the expression to zero, it allows for any value of $c$, and thus the nullclines for nonzero c are horizontal lines in the $c-m$ plane.
    The fundamnetal theorem of algebra promises us at least one real root for the cubic. For three distinct real roots, the discriminant of the cubic must be positive. The discriminant of a cubic polynomial $ax^3 + bx^2 + c x + d$ is of the form:
\begin{align} \Delta &= b^2 c^2 - 4 a c^3 - 4 b^3 d - 27 a^2 d^2 + 18 a b c d \\[0.5em] & \hspace{2em} a = \alpha_c, \hspace{0.5em} b = \gamma-\beta_0, \hspace{0.5em} c = \alpha_c, \hspace{0.5em} d = \gamma\\[0.5em] \Delta &= (\gamma - \beta_0)^2 \alpha_c^2 - 4 a_c^4 - 4(\gamma - \beta_0)^3\gamma - 27 a_c^2 \gamma^2 + 18 \alpha_c^2 (\gamma-\beta_0)\gamma \end{align}

When $\Delta > 0$: \begin{align} \alpha_c^2 [(\gamma - \beta_0)^2 - 27 \gamma^2 18 \gamma (\gamma - \beta_0)] - 4 \gamma (\gamma - \beta_0)^3 - 4 \alpha_c^4 &> 0 \\[0.5em] \alpha_c^2 [(\gamma-\beta_0)(\gamma-\beta_0+18\gamma)-27\gamma^2] - 4 \gamma (\gamma - \beta_0)^3 - 4 \alpha_c^4 &> 0 \\[0.5em] \alpha_c^2 [-8\gamma^2 + \beta_0^2 - 20 \gamma\beta_0] - 4 \gamma (\gamma - \beta_0)^3 - 4 \alpha_c^4 &> 0 \\[0.5em] \boxed{\alpha_c^2 [\beta_0^2 - 20 \gamma\beta_0 - 8\gamma^2] + 4 \gamma (\beta_0 - \gamma)^3 - 4 \alpha_c^4 > 0}\\[0.5em] \end{align} I believe there is a typo in the problem statement, (the last term $4\alpha_c^4$ should have a power of 4 and not 2). I could be wrong though (I pulled the determinant equation off of Wikipedia), but I suspect it's a typo otherwise the squared term would be grouped together and simplified.

c) Nullcline for $\cfrac{\mathrm{d}m}{\mathrm{d}t} = 0$

Rewriting \begin{align} \cfrac{\mathrm{d}m}{\mathrm{d}t} = 0 &= I_0 + c - \alpha_0 c \cfrac{m^2}{1+m^2} - m \\[0.5em] &0 = I_0(1+m^2) + c(1+m^2) - \alpha_0 c m^2 - m(1+m^2)\\[0.5em] &0 = -m^3 + I_0 m^2 - m + I_0 + c(1 + m^2 - \alpha_0 m^2)\\[0.5em] &c = \cfrac{m^3 - I_0 m^2 + m - I_0}{1+m^2(1-\alpha_0)}\\[0.5em] &\boxed{c = \cfrac{(m-I_0)(1+m^2)}{1+m^2(1-\alpha_0)}}\\[0.5em] \end{align}

This function is undefined when the denominator approaches zero, so: $$m_{st}^2 \neq \cfrac{1}{\alpha_0-1}$$ Furthermore, for c > 0 requires either that both numerator and denominator are positive, or both are negative. Note that when $\alpha_0 \leq 1$, the denominator is always positive. So we can break up our cases like so:

$\alpha_0 \leq 1$:

When $\alpha_0$ is less than 1, there is a sharp lower bound for steady state $m$. $$ \boxed{m_{st} \geq I_o }$$

$\alpha_0 > 1$:

When $\alpha_0$is greater than 1, either the first inequality holds or the second. The steady state concentration of the cytokine is sandwiched between these two parameters. \begin{align} \boxed{ \hspace{3em} I_0 < m_{st} < \sqrt{\cfrac{1}{\alpha_0-1}} \\[0.5em] \sqrt{\cfrac{1}{\alpha_0-1}} < m_{st} < I_0 } \end{align}

Let's rewrite our original ODE and collect the coefficient of $c$ for some insight: $\cfrac{\mathrm{d}m}{\mathrm{d}t} = I_o + c\left(1-\alpha_0\cfrac{m^2}{1+m^2}\right) - m $.
It now becomes more clear what the threshold at $\alpha_0 = 1$ means; $m$'s response to concentrations of $c$ is either capable of being a net positive or negative (the hill function activation term is bounded by one and reaches one for m sufficiently large). So when $\alpha_0 > 1$, dm/dt is cleared (taken in by cells) faster than the per cell rate at which it is being secreted (the original dimensionalized inequality looks at $\alpha_0$ vs. $\beta_1$).

Thus, if $\alpha_0$ is high, cytokine uptake is so high it sets what the upper bound of what $m_{st}$ could ever be when $I_0$, the basal level of production, is small. Conversely, when $I_0$ is relatively high, $m$ can recover from uptake by the cells and is relatively unaffected. But because the cells are also secreting cytokine, there is an inherent lower bound within $\alpha_0$.

Analytically, we can think of $m_{st} (\alpha_c, \gamma, \beta_0)$, and $c_{st}(m_{st}, I_0, \alpha_0)$, where both steady states are determined purely by the adjacent parameters. From a design standpoint, if we wanted two stable states we want to have an $I_o$ and $\alpha_0$ that satisfies the relation above.

d) Fixed points

  1. When c = 0, $\mathrm{d}c/\mathrm{d}t = 0$ is already satisfied, and $dm/dt=0$ when the cubic $m^3 - I_0 m^2 + m - I_0 = 0$, or with $C$ and $D_0$ defined above,
    \begin{align} m = \cfrac{1}{3}\left(I_0 - C - \cfrac{D_0}{C}\right) \end{align}
  2. Other two fixed points can be found by first looking at $dc/dt = 0$ and setting $\alpha_c m_{st}^3 + (\gamma-\beta_0)m_{st}^2 + \alpha_c m_{st} + \gamma = 0$. Then, finding $c_{st}$ from $c_{st} = (m^3 - I_0 m^2 + m - I_0)/(1+m^2(1-\alpha_0))$.
    Use parameters $\beta_0 = 9.78, \gamma = 2.14, I_0 = 24.67, \alpha_c = 1.06, \alpha_0 = 11.21$.
  3. Plotting nullclines...

The log toggle does not add much, it just helps zoom in on the bottom right fixed point. (note that $c = 0 \rightarrow \log c = -\infty$, so the logged plot only shows one stable fixed point and is limited in this sense)

  1. Magnitudes of $I_0$ and $\alpha_0$:
    I am not sure why this question is placed here. I also did a little more exploration in part (h) later in the notebook, and so there will be more coming then, but for now, I've copied and pasted my analysis from part (c) here:

    Let's rewrite our original ODE and collect the coefficient of $c$ for some insight: $\cfrac{\mathrm{d}m}{\mathrm{d}t} = I_o + c\left(1-\alpha_0\cfrac{m^2}{1+m^2}\right) - m $.
    It now becomes more clear what the threshold at $\alpha_0 = 1$ means; $m$'s response to concentrations of $c$ is either capable of being a net positive or negative (the hill function activation term is bounded by one and reaches one for m sufficiently large). So when $\alpha_0 > 1$, dm/dt is cleared (taken in by cells) faster than the per cell rate at which it is being secreted (the original dimensionalized inequality looks at $\alpha_0$ vs. $\beta_1$).
    Thus, if $\alpha_0$ is high, cytokine uptake is so high it sets what the upper bound of what $m_{st}$ could ever be when $I_0$, the basal level of production, is small. Conversely, when $I_0$ is relatively high, $m$ can recover from uptake by the cells and is relatively unaffected. But because the cells are also secreting cytokine, there is an inherent lower bound within $\alpha_0$. From a design standpoint Analytically, we can think of $m_{st} (\alpha_c, \gamma, \beta_0)$, and $c_{st}(m_{st}, I_0, \alpha_0)$, where both steady states are determined purely by the adjacent parameters. From a design standpoint, if we wanted two stable states we want to have an $I_o$ and $\alpha_0$ that satisfies the relation above.

e) separatrix & f) phase portrait

Linear stability matrix \begin{align} A = \begin{pmatrix} \beta_0 \cfrac{m^2}{1+m^2} - \gamma - \alpha_c m, & \beta_0 \cfrac{2cm}{(1+m^2)^2} - \alpha_c c\\ 1 - \alpha_0 \cfrac{m^2}{1+m^2}, & - \alpha_0 \cfrac{2cm}{(1+m^2)^2} - 1 \end{pmatrix} \end{align}

Since this code is a direct copy-paste of what's above, except the phase portrait which is really just calling biocircuits, I've abstracted the code away into a module called utilities.py. Let's call the function and build our plot! (Note that to view the code of the function in notebook can run the following line: utilities.part_ef_plotter??)

g) Temporal profiles above and below the separatrix

So we see above the separatrix, we converge to the concentrations (c, m) of the stable fixed point (0, 24.7), and below the separatrix, we converge to the point (10.8, 0.64). Looking at the shape of the separatrix, we see that to reach the homeostatic stable state (so that cytosine concentration doesn't just collapse), the initial m generally has to be below ~2-3x the initial c.
The large divisive basins provide a threshold for what we colloquially refer to as homeostasis, where perturbations within the basin below the separatrix demonstrates a robustness to the initial starting concentrations, only sensitive along the separatrix itself (as demonstrated in the plot above).
Otherwise the system's steady state remains robust to a wide range of initial conditions (shown below). Here, we take our separatrix, then shift sep_m by 1 below the separatrix line, and plot all the traces of this trajectory. Since the code is really just the same as what's shown above, the function is stored in utilities.py as part_g_plotter().

h)

I was going to try a stochastic simulation next, but I got distracted and found this lil easter egg in the comments of the problem statement!

We are cutting this part of the problem just to keep the homework from getting too long.

h) Play with parameter values and investigate conditions for bistability. You already have much the machinery in place to do this, since you can already computed the fixed points. For bistability, we must have three total fixed points (two in addition to the $C = 0$ fixed point that always exists). You can use this as a criterion to find ranges of parameter values for which bistability exists. With respect to which parameter values is bistability most robust? How about least robust?

I abstracted a lot of the code away into a module utilities.py. The major addition is defining a function draw_fp() that will responsively plot different glyphs for different stabilities, defining a dashboard-creator param_plotter() that will plot our nullclines as we toggle our parameters. Here is the result:

Note that as soon as roots become complex, the points disappear (try-except-pass to the rescue!) MORE ANALYSIS

where we are...

So we make the following conclusion from our mathematical analysis in parts (a) - (g) and the visualizations above. The existence of $m_{st}$ depends on the existence of positive real roots in the cubic that is parametrized by three parameters: $(\alpha_c, \gamma, \beta_0)$. We saw in part (b) that a necessary condition is that $\gamma < \beta_0$ and that the discriminant, a quadratic in $\alpha_c^2$, be positive. Once the root of $m_{st}$ is found, $c_{st}$ can be found by plugging in $m_{st}$ to an equation parametrized by $(\alpha_0, I_0)$. The positive nature of $c_{st}$ is then dependent on the inequalities we set up in part (d.iv).

where we're going...

The authors of the paper provided values that we used for parts (a)-(g). Now, we will broaden our search:

Even with the json files, the slider is really slow on my computer (150x150 grid for 2 plots with a slider that takes on 20 values means 900,000 points in total) so I tried turning it into a movie—hopefully it works!

Here is the dashboard:

So we see that we do not want $\gamma$ to be too high, otherwise we won't be able to find an $I_0$ or $\alpha_0$ at a given $\alpha_c$, $\beta_0$. But once we satisfy that condition, we can have a stable fixed point when both $I_0$ and $\alpha_0$ are of similar orders of magnitude. The block-like shape emerges given the chained inequalities we set up earlier, where our root is satisfied either with $I_0$ as a lower bound (upper right block), or as an upper bound (lower left block), and vice veras for $\alpha_0$.

To answer the original question posed in the prompt, the system appears most robust to $Io, \alpha_0$, given the broad ranges they can take on. In regions of small $\alpha_c$, $\beta_0$ can take a range of values (big red block at $\gamma = 5.5$). I am hesitant to say that the system is most sensitive to $\gamma$, since that is the parameter I have chosen as my slider; it could easily be the case we lose stability in the right plot when looking along other dimensions. However, it is certainly the relation between $\beta_0$ and $\gamma$ that results in greater sensitvity, for we see that $alpha_c$ can take on values from $10^{-3} - 10^2$, whereas $\beta$ and $\gamma$ are both sandwiched for the most part between $10^1 - 10^2$. (note that although the left plot is a square, the scale of the axes are not, so the y-chunk is narrower than the x-chunk in actual parameter space).

Thank you much for reading!! This is all for 7.1.

Note: Is it ok that I am putting functions in a separate module? I just find the scrolling to be so overwhelming that I want the focus to be on my analysis and the visuals instead of the enormous amount of semi-repetitive code. Please let me know if this is counterproductive, as I can understand why it is important to be able to look at and interact with the code to see if it is functioning as it claims. That being said, I've included an "Appendix" for my functions in utilities, and the code that was used to create the .json files

Appendix A)

Appendix B)

Here is the code that was used to create the abg.json and aig.json files.

Appendix C)

~ Poem? ~

Rise and shine, my Cytokine.
Wine and dine,
down my spine,

align your signs—
for bistability is divine!